Loading Data
epi_rpca <- readRDS(file.path(PATH, "data/sc/epi_rpca_subset.rds"))
epi_celltypes <- readRDS(file.path(PATH, "brca_atlas_validation/data/sigs/epi_markers.rds"))
epi_barkley <- readRDS(file.path(PATH, "data/signatures/cancer_epithelial/barkley_natgen_2023/barkley_sigs.rds"))
epi_gavish <- readRDS(file.path(PATH, "data/signatures/cancer_epithelial/gavish_nature_2023/gavish_sigs.rds"))
epi_xu <- readRDS(file.path(PATH, "data/signatures/cancer_epithelial/xu_cellreports_2024/xu_sigs.rds"))
epi_celltypes <- lapply(epi_celltypes, function(x) x[x %in% rownames(epi_rpca)])
epi_barkley <- lapply(epi_barkley, function(x) x[x %in% rownames(epi_rpca)])
epi_gavish <- lapply(epi_gavish, function(x) x[x %in% rownames(epi_rpca)])
epi_xu <- lapply(epi_xu, function(x) x[x %in% rownames(epi_rpca)])
epi_all <- c(epi_celltypes, epi_barkley, epi_gavish, epi_xu)
Annotations
epi_rpca$celltypist_epi_only <- with(epi_rpca@meta.data,
case_when(celltypist_broad == "Epithelial" ~ celltypist_pred,
.default = celltypist_broad))
epi_rpca$singler_epi_only <- with(epi_rpca@meta.data,
case_when(singleR_broad == "Epithelial" ~ singler_pred,
.default = singleR_broad))
epi_rpca$author_epi_only <- with(epi_rpca@meta.data,
case_when(author_broad == "Epithelial" ~ celltype_new,
.default = author_broad))
p1 <- DimPlot_scCustom(epi_rpca,
reduction = "umap.rpca",
group.by = "RNA_snn_res.0.2",
label = TRUE,
label.size = 3,
raster = FALSE,
repel = TRUE) + labs(title = "Clustering (0.2)", x = "UMAP 1", y = "UMAP 2")
p2 <- DimPlot_scCustom(seurat_object = epi_rpca,
group.by = "celltypist_epi_only",
reduction = "umap.rpca",
label = TRUE,
label.size = 1.5,
label.box = TRUE,
repel = TRUE,
raster = FALSE) +
labs(title = "Celltypist", x = "UMAP 1", y = "UMAP 2") +
theme(legend.position = "none")
p3 <- DimPlot_scCustom(seurat_object = epi_rpca,
group.by = "singler_epi_only",
reduction = "umap.rpca",
label = TRUE,
label.size = 1.5,
label.box = TRUE,
repel = TRUE,
raster = FALSE) +
labs(title = "SingleR", x = "UMAP 1", y = "UMAP 2") +
theme(legend.position = "none")
p4 <- DimPlot_scCustom(epi_rpca,
reduction = "umap.rpca",
group.by = "author_epi_only",
label = TRUE,
label.size = 1.5,
label.box = TRUE,
repel = TRUE,
raster = FALSE) +
labs(title = "Author", x = "UMAP 1", y = "UMAP 2") +
theme(legend.position = "none")
combined <- cowplot::plot_grid(p1, p2, p3, p4, ncol = 2, nrow = 2)
ggsave(plot = combined, filename = file.path(PATH, "results/umaps/epi_rpca_annot_02.png"), height = 8, width = 9)
combined
Hormone/HER2 Expression
p1 <- DimPlot_scCustom(epi_rpca,
reduction = "umap.rpca",
group.by = "RNA_snn_res.0.2",
label = TRUE,
label.size = 3,
raster = FALSE,
repel = TRUE) + labs(title = "Clustering (0.2)", x = "UMAP 1", y = "UMAP 2")
p2 <- DimPlot_scCustom(epi_rpca,
reduction = "umap.rpca",
group.by = "subtype_new",
label = TRUE,
label.size = 3,
raster = FALSE,
repel = TRUE) + labs(title = "Subtype", x = "UMAP 1", y = "UMAP 2")
p3 <- FeaturePlot_scCustom(epi_rpca,
reduction = "umap.rpca",
features = "ESR1")
##
## NOTE: FeaturePlot_scCustom uses a specified `na_cutoff` when plotting to
## color cells with no expression as background color separate from color scale.
## Please ensure `na_cutoff` value is appropriate for feature being plotted.
## Default setting is appropriate for use when plotting from 'RNA' assay.
## When `na_cutoff` not appropriate (e.g., module scores) set to NULL to
## plot all cells in gradient color palette.
##
## -----This message will be shown once per session.-----
p4 <- FeaturePlot_scCustom(epi_rpca,
reduction = "umap.rpca",
features = "ESR2")
p5 <- FeaturePlot_scCustom(epi_rpca,
reduction = "umap.rpca",
features = "PGR")
p6 <- FeaturePlot_scCustom(epi_rpca,
reduction = "umap.rpca",
features = "ERBB2")
combined <- cowplot::plot_grid(p1, p2, p3, p4, p5, p6, ncol = 3, nrow = 2)
ggsave(plot = combined, filename = file.path(PATH, "results/umaps/epi_rpca_subtype_02.png"), height = 8, width = 12)
combined
AUC scores
if(do_save){
K2_sce <- as.SingleCellExperiment(epi_rpca)
cells_auc <- AUCell_run(K2_sce, epi_all)
saveRDS(cells_auc, file.path(PATH, "data/auc/epi_auc.rds"))
} else {
cells_auc <- readRDS(file.path(PATH, "data/auc/epi_auc.rds"))
}
Dot plots
all.equal(epi_rpca %>% colnames, cells_auc@assays@data$AUC %>% colnames)
## [1] TRUE
auc_assay <- CreateAssay5Object(counts = cells_auc@assays@data$AUC)
## Warning: Data is of class matrix. Coercing to dgCMatrix.
auc_seurat <- CreateSeuratObject(auc_assay,
meta.data = epi_rpca@meta.data)
Idents(auc_seurat) <- auc_seurat$RNA_snn_res.0.2
Clustered_DotPlot(seurat_object = auc_seurat, features = names(epi_all), k = 19)
## Warning: Feature list contains duplicates, making unique.
## Warning: No layers found matching search pattern provided
## Warning in FetchData.Assay5(object = object[[DefaultAssay(object = object)]], :
## data layer is not found and counts layer is used
## [[1]]
##
## [[2]]